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ABSTRACT 

We have performed three-dimensional numerical simulations of accretion discs in a 
close binary system using the Smoothed Particle Hydrodynamics (SPH) method. Our 
results show that, contrary to previous claims, 3D discs do exist even when the specific 
heat ratio of the gas is as large as 7 = 1.2. Although the disc is clearly more spread 
in the z-direction in this case than it is for the quasi-isothermal one, the disc height 
is compatible with the hydrostatic balance equation. 

Our numerical simulations with 7 = 1.2 also demonstrate that spiral shocks exist 
in 3D discs. These results therefore confirm previous 2D simulations. 
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1 INTRODUCTION 

Numerical studies of accretion discs have been mostly re- 
stricted to 2-dimensional cases, due to computing time lim- 
itations. Among many things, these 2D simulations have 
shown that spiral shocks appear in inviscid discs (Sawada et 
al. 1987, Rozyczka & Spruit 1989; Spruit 1989). These spiral 
shocks have been invoked as an alternative to the a- viscosity 
theory for transfering angular momentum outward, thereby 
allowing accretion. It is however still unclear whether these 
spiral shocks are effective enough to explain observations. 

Recently, some 3D simulations have been carried out. 
Except for Sawada & Matsuda (1992), they were all per- 
formed using particle methods (Molteni, Belvedere & Lan- 
zafame 1991; Hirose, Osaki & Mineshige 1991; Lanzafame, 
Belvedere & Molteni 1992, 1993; Meglicki, Wickramasinghe 
& Bicknell 1993; Whitehurst 1994; Simpson 1995; Armitage 
& Livio 1996). 

Sawada and Matsuda (1992) used the Eulerian grid 
based Roe upwind TVD scheme to calculate the case 7 = 1.2 
and mass ratio M1/M2 = 1.0 up to half an orbital period. 
They found the accretion disc to be not axisymmetric but 
characterized by a pair of spiral shocks, although their com- 
putational time was not long enough to reach a steady state. 

Hirose et al. (1991) used the "sticky" particle method of 
Lin & Pringle (1976). In this method, the effect of pressure 
is neglected but not viscous interactions, and the energy 
generated is assumed to be radiated away instantaneously. 
Hirose et al. (1991) considered two values of the mass ratio, 
0.15 and 1. They observed that the vertical height of the 
disc, H , is generally much larger than the value expected 
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from the hydrostatic balance equation. The height of the 
disc also appeared to be a variable function of the orbital 
phase and of the distance (r). They attributed this to the 
interaction between an incoming stream and the gas in the 
disc. The ratio H/r they obtained was typically between 10 
and 15 %. 

Molteni et al. (1991) carried out 3-D simulations of the 
formation and the evolution of polytropic accretion discs 
using the SPH method. They made calculations for three 
values of the polytropic index, that is 7 = 1.01, 1.1 and 
1.2, using 9899, 2372 and 1579 SPH particles respectively. 
They concluded that the formation of a disc is inhibited 
for 7 > 1.2 and that the thin disc approximation is con- 
firmed only for low 7 values. Later, Lanzafame et al. (1992) 
performed calculations with a mass ratio M1/M2 = 1.625. 
They concluded that disc formation is inhibited for 7 > 1.1, 
a value still lower than that found in the previous paper. 

Meglicki et al. (1993), in an attempt to study large- 
scale eddies, also used SPH but with a very low value for 
the artificial viscosity. Their model had been evolved only 
for 1.7 orbital periods but showed the presence of turbulence 
which may play a role in transporting angular momentum. 
Turbulent eddies were indeed seen up to scales approaching 
half the thickness of the disc. 

Whithehurst (1994) studied the eccentric disc model of 
superhump formation in SU UMa stars using a free Lagrange 
method which approximates the gas as a collection of fixed- 
mass cells with variable volumes. He typically used about 
60 000 cells in his 3D calculation of accretion disc with a 
mass ratio smaller than 0.333 . An ad-hoc artificial viscosity 
is included in his method to allow for angular momentum 
transfer in a way comparable to the standard a-disc the- 
ory. Energy dissipated was allowed to radiate away instan- 
taneously. 
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Simpson (1995) calculated the case 7 = 1.01 and 
M1/M2 = 12.5 until about 50 orbital periods with SPH, 
using 5 000 particles. No eccentric modes were seen to de- 
velop. However, spiral density waves were visible for most of 
the later orbits. 

Armitage and Livio (1996) investigated the collision be- 
tween the stream from the inner Lagrange point, Li, and the 
accretion disc using the SPH method, in a binary system 
with a mass-ratio of 0.2. They used an isothermal equation 
of state and set-up initially the disc as an annulus of 30,000 
particles. They observed that significant amounts of stream 
material overflowed or bounced off the edge of the disc and 
flowed over the disc surface toward smaller radii. This, they 
believe, is a likely cause of the absorption dips observed in 
some nearly edge-on low-mass X-ray binaries. 

One of the most striking facts from most of these parti- 
cle calculations is their apparent inability to generate spiral 
shocks in the accretion discs. This may bo related to the 
fact that, except for Molteni et al. (1991) and Lanzafame 
et al. (1992), the authors considered isothermal or pseudo- 
isothermal cases or they used a large amount of viscosity in 
their method. Therefore, these calculations do not constrain 
very much the formation of spiral shocks in 3D. On the other 
hand if spiral shocks were to exist in 3D, the calculations of 
Molteni et al. (1991) and Lanzafame et al. (1992) should 
have shown them. However, as already discussed in Boffin, 
Yukawa & Matsuda (1997; hereafter Paper I), their use of a 
very small number of particles in the 7 — 1.1 and 1.2 cases, 
of an arithmetic mean for the pressure gradient as well as of 
a constant smoothing length, may lead to incorrect results. 

We would like therefore, in this paper, to address two 
questions : (i) is disc formation really inhibited for value of 
7 > 1.1 ? (ii) Can wo extend the results of Paper I and also 
observe spiral shocks in 3D ? 



2 METHOD 

We use the SPH method (see e.g. Monaghan 1992 for a re- 
view) , which is rather well suited to the study of mass trans- 
fer in a binary system. The details of our numerical method 
are described in Paper L 

We choose the following units: orbital separation A, 
total mass M = Mi -|- M2 and orbital velocity Vorb = 
[G{Mi+M2)/AY^^. The orbital period is therefore normal- 
ized to 27r. We will study only cases with mass ratio equal 
to unity, i.e. Mi = M2 = 0.5. 

We assume a non- viscous gas and use the standard SPH 
artificial viscosity with a — 1.0 and /3 — 2.0. We do not take 
into account any radiative cooling effect. We work in the 
non-inertial frame of reference where both stars are at rest. 
The primary is centered at the origin, the secondary at (- 
1,0,0). 

We carry out two kinds of simulations : with and with- 
out mass supply from the secondary star. In both cases, the 
computing region was a sphere of 0.55 around the accreting 
star whose radius is 0.03. 

In the Roche lobe overflow(RLOF) case, gas from the 
secondary star flows with the velocity of sound through the 
Li point toward the primary. The sound speed of the gas 
is 0.1 and its density is 1. SPH particles at the Li point 
are ordered on a square lattice of width 0.04. Throughout 
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Figure 1. Evolution of the disc mass as a function of time, for 
the 3 simulations of Sect. 3.1 



the calculation, we keep a constant value of the injection 
rate, 1.6 x 10^ . The initial value of the smoothing length 
is 0.0026 (7 = 1.2), 0.00308 (7 = 1.1) or 0.005 (7 = 1.01). 
At the end of the simulation, it ranges from 0.0026 to 0.85, 
with a mean value of 0.0067 (7 =1.01) or 0.0094 (7 =1.2). 

In the second kind of simulation, we calculate the evo- 
lution of the disc from an initial state in which each parti- 
cle has been given a Keplerian velocity. Here, only the case 
7 = 1.2 has been investigated. The radius of the initial disc 
is 0.30 and the thickness is 0.08. The gas in the disc has 
also a sound speed of 0.1 and a density of 1. The vertical 
structure of the disc is obtained by placing 9 layers - each 
separated by 0.01 - on top of each other. This ensures that 
the disc is resolved vertically. Very quickly, the disc adjusts 
itself so that the particles distribution in the 2-direction is 
gaussian (see also below). Each layer consists of an annu- 
lus between r = 0.05 and r = 0.30. Particles are placed on 
concentric rings such that the mean distance between par- 
ticles is 0.01. Initially, the smoothing length is 0.011 and 
there are 31 000 particles. At the end of the simulation, at 
about T=2.3, the moan value of the smoothing length for 
the remaining 18,000 particles is 0.026 . 



3 RESULTS 

3.1 Roche lobe overflow 

Wo have calculated three cases corresponding to 7 = 1.2,1.1 
and 1.01, until T ~ 14 (that is, a little more than 2 or- 
bital periods), T ~ 19 and T ~ 39, respectively. As can 
be seen from Fig. 1, which shows that the disc masses are 
still increasing, we have not yet reached a steady state. It 
is interesting to note that the prcssurclcss calculations of 
Hirose et al. (1991) reached a steady state after only about 
100 orbital periods, while Molteni et al. (1991) get a station- 
ary state after three orbital periods or less. Our 2D results 
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Figure 2. The time evolution of the particle location map in the orbital plane in the case of 7 = 1.2. 



Figure 2. Continued 




Figure 3. Gray scale plot of density in the x — y plane at time 
T ~ 14 in the case of 7 = 1.2. 



Figure 4. The location of particles in the x — z plane at time 
T ~ 14 in the case of 7 = 1.2 ; 13730 particles out of a total of 
61369 are shown. 



(see Paper I) show also that the mass of the disc is station- 
ary only after several tens of orbital periods. However, from 
these 2D results, as well as from the fact that the gross fea- 
tures in the flow do not change very much, we are confident 
that we have almost obtained the final structure of the disc. 

The time evolution of the flow is shown for the 7 = 1.2 
case in Fig. ^. Only particles within 0.02 of the orbital plane 
are represented. It is clear that the structure of the flow 
does not change after, say, one orbital period. This was also 
observed in our 2D results. But, although the basic pattern 
is unchanged, the density in the disc increases continuously 
as a result of an increase of the number of particles. From 
Fig. |2| it is also obvious that spiral shocks, similar to those 
seen in 2D, are present in the orbital plane. 

We show our final results for the 7 — 1.2 case in Fig. ^ 
and following. At this stage, we have 61369 particles. The 
radius of the disc is approximately 0.2, similar to the 2D 
results. Fig. ^ clearly shows that the disc is asymmetric in 
the orbital plane: the upper half is more spread out than the 
lower one. The disc is however totally symmetric perpendic- 
ular to the equatorial plane, as shown by Fig. ^. 

In order to help understand the fiow structure, we anal- 
yse now in more detail the velocity profile. First, in Fig. ^ 
the azimuthal and radial components of the velocity are 
shown as a function of the azimuthal angle, 9, for three dif- 
ferent ranges in radial distance. For the azimuthal velocity, 
the two horizontal solid lines represent the Keplerian ve- 
locity corresponding to the inner and outer radius of each 
region. In the radial velocity plots, the dashed lines indicate 



the approximate positions of the shock. Figure ^ clearly 
demonstrates the effect of the shock : the azimuthal com- 
ponent of the velocity decreases while the radial component 
changes sign. In front of the shocks it has a plus sign, mean- 
ing gas is accelerated outward, while behind them it has a 
minus sign, implying an inward acceleration of the gas. It 
must be noted that in the 0.05 < r < 0.10 range, the veloc- 
ity profile is somewhat less clear. This is due to the inner 
vacuum boundary condition which reduces the number of 
particles present in this region. 

A definitive confirmation of the presence of shocks could 
be done by checking that the normal component of the ve- 
locity with respect to the shock shows a transition from su- 
personic to subsonic flow. This, however, requires the knowl- 
edge of the pitch angle of the spiral shocks, which is difficult 
to estimate. We have thus to limit ourselves to analyse the 
azimuthal and radial Mach numbers. The mean azimuthal 
Mach number is around 5 which, if the hydrostatic balance 
is valid, would give a disc height to disc radius ratio of 0.2. 
This is indeed the case as shown below. On the other hand, 
the radial Mach number is - except for a few particles - al- 
ways smaller than 1. It is however interesting to note that 
the shock clearly reduces the radial component of the veloc- 
ity, bringing it from near supersonic value to zero. 

Figure ^ is a gray scale plot of the radial Mach number. 
The presence of the two-armed spiral clearly demonstrates 
the effect of the shock on the velocity. It can also been seen 
that the radial Mach number changes sign along the spiral 
shock surfaces. 

We have also estimated the azimuthally averaged angu- 
lar advected momentum fiux as well as the torque density 
between the companion star and the disc at our final stage 
for the case 7 = 1.2. This is shown in Fig. |^. 

We can see, in the 0.06 < r < 0.15 range, that the 
torque is negative, meaning that the disc loses its own angu- 
lar momentum and gives it to the orbital angular momentum 
of the binary system. For r < 0.06, the effect of the inner 
boundary (vacuum) condition is seen, namely that the gas is 
accelerated inward. As the radius of the disc is roughly 0.2 
(see Fig. 3), the region with r > 0.2 is not very meaningful 
to consider. However, the same tendency as in the disc re- 
gion is seen, i.e. the negative torque density contributes to 
the angular momentum transfer. 

In Fig. 1^ and ^ we show our final results for the 7 = 
1.1 and 1.01 cases, respectively. In both cases, it is more 
difficult to see any spiral shock although the 7 = 1.1 figure 
shows a build-up of particles in the disc. These results are 
rather similar to those we obtain in 2D (Paper I) . When we 
consider the mean density as a function of the distance to 
the primary, the three cases show a similar pattern : density 
increases sharply until about r ~ 0.1, then decreases rather 
quickly. In the case of 7 = 1.01, at r = 0.2, the mean density 
is only 10 % of the maximum value, while it is about 40 % for 
7 — 1.2. The maximum value of the mean density itself is, 
however, very dependent on 7 : 1.8, 0.4 and 0.09 for 7 = 1.01, 
1.1 and 1.2, respectively. 
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From Fig. we can have an idea of the disc height 
in the 7 — 1.2 case. Although a relatively large spread 
can be seen, the disc appears rather well bound. We thus 
believe that even in the case of a large polytropic index, 
a disc forms. Using the verified assumption that the par- 
ticles' distribution in the z-direction is a Gaussian, i.e. 
dn{z) - exp{-z^/2H'^)dz, we estimated the mean height 
of the disc as a function of the distance in our 3 simulations. 
A linear trend was found. The ratio H/r is about 0.09, 0.12 
to 0.2 and 0.2 to 0.3 for 7 = 1.01, 1.1 and 1.2, respectively. 
Thus, the disc height is about 2 to 3 times larger for 7 = 1.2 
than it is for 7 = 1.01. However, one should note that in 
all cases, H almost exactly equals the value expected from 
the hydrostatic balance, Cs/Q.. Therefore, our discs appear 
to be stable. Of course, inclusion of radiation loss would be 
required to simulate real discs. This would reduce the sound 
speed close to the star and hence make the disc thinner. 
Concerning the disc height, we should also mention that we 
could not observe any variation with the orbital phase as, 
for example, Hirose et al. (1991) did. Nor did we observe 
asymmetry of the flow for particles outside the orbital plane 
as reported by Molteni et al. (1991) and Lanzafame et al. 
(1992). 

Finally, we turn to the accretion rate, shown in Fig. |lo[ 
As already stressed at the beginning of this section, a steady- 
state has not been reached. This can also be seen in Fig. |lO[ 
For now, we observe that the mass accretion rate is slightly 
dependent of the polytropic index : the smaller 7 is, the 
smaller the accretion rate is. Roughly, the ratio between the 
injection rate and the accretion rate is on the order of 50 % 
to 60 %. Paper I has however shown that the accretion rate is 
a function of the numerical resolution. We cannot therefore 
claim to provide the definitive number. Also from Paper 
1, it is not clear whether the mass accretion rate reported 
here is entirely due to the presence of the spiral shocks as a 
large fraction may be due to the numerical viscosity. Indeed, 
in Paper I, we investigated a model in which we did not 
take into account the presence of a companion, hence we 
studied an isolated accretion disc. In this case, accretion 
is only due to viscous effect as no tidal torque is present. 
However, we noted that the accretion rate so obtained was 
of the same order of magnitude as when a companion is 
present. It is thus difficult to disentangle the effect due to 
numerical viscosity and to the companion. Another problem 
arises from our inner boundary condition. Here, we simply 
remove particles when they enter the boundary. The vacuum 
so created results in a spurious force which increases the 
accretion rate. We therefore suspect that the difficulty of 
obtaining a quantitative estimate of the accretion rate is 
related to two effects : the SPH articifial viscosity term and 
the inner boundary condition. We are presently trying to 
solve these problems for further investigations. 

It may, however, be interesting to note that the ratio 
of 50-60 % is rather similar to the value obtained in our 2D 
results. On the other hand, the escape rate is much smaller 
than the value obtained in the 2D cases at the same time. 
This can be easily understood : as the gas can now expand 
in the z-direction, it is more difficult for it to fill the Roche 
lobe and therefore escape. 

As we have not yet reached a steady state, this may 
only be a temporary effect. One should also note that, as in 
the 2D results, the escape rate is a function of 7. 
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Figure 6. Gray scale plot of the radial Mach number in the x — y 
plane for the 7 = 1.2 case at time T ^ 14. 
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Figure 7. The azimuthal averaged advected angular momentum 
flux(a) and tidal torque density(b) as a function of radius for the 
case 7 = 1.2 at time T ~ 14. The average was done on particles 
having \z\ < 0.02 only. 



Figure 8. The particle location map in the orbital plane in the 
case 7 = 1.1 at time T ~ 19 : 34779 particles out of a total of 
55370 are shown. 



Figure 9. The particle location map in the orbital plane in the 
case of 7 = 1.01 at time T ~ 39 : 23960 particles out of a total of 
26375 are shown. 
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Figure 5. Azimuthal plots of the azimuthal and radial components of the velocity for the case 7 = 1.2 at time T ^ 14. Only particles 
with \z\ < 0.02 are shown. The two solid lines in the azimuthal velocity plots represent the Keplerian velocity corresponding to the inner 
and outer radius of each region. The dashed lines indicate the approximate positions of the shock. 



3.2 Keplerian disc 

Figure ^ shows the particle distribution in the orbital plane 
in the case of no mass inflow. Although the simulation could 
be followed for only a short time, the spiral shock is again 
clearly seen, as well as the elongation of the disc in the y- 
direction. From a close inspection of the particle distribu- 
tion, it appears that the spiral pattern does not change in 
the 2-direction. We also estimated the disc height using the 
verified assumption of a Gaussian distribution. The mean 
disc height, H , is phase independent and scales almost lin- 
early with the distance to the primary star in the equatorial 
plane, r, in the following way : 

= 0.13 r 0.011. (1) 



Although this shows that H may take values between 
about 20 and 35 % of r, it almost exactly equals the value 
expected from the hydrostatic balance. Thus, in this case 
also, the disc height is only due to the pressure. Pressure is 
also responsible for the spread in velocity around the Kep- 
lerian value. We can see no difference in the disc height or 
in the velocity distributions of this simulation and the one 
described in the previous section. We can therefore not con- 
firm the claim of Hirose et al. (1991) that the interaction 
between the Li flow and the disc is responsible for a large 
and spatially-variable disc height. We recall that Hirose et 
al. (1991) did not include pressure in their calculations. 
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Figure 10. Time history of the mass accretion rate for 7 = 1.01 
(full curve), 1.1 (broken curve) and 1.2 (dotted curve). The unit 
is the mass injection rate, 0.00016. 



SPH method (which includes the eflFect of pressure), with 
about (5 — 6) X 10* particles and a variable smoothing length. 

(2) In the 7 = 1.2 case, spiral shocks are seen. These 
spiral shocks have a clear effect on the velocity of the gas. 
In the 7 = 1.1 and 7 = 1.01 cases, although spiral shocks 
are not clearly seen, a build-up of density is present. These 
results axe in agreement with our 2D results presented in 
Paper I. This gives strong support to the accuracy of 2D 
simulations. 

(3) Spiral shocks wore also observed in discs without 
mass supply from the inner Lagrange point. In this case, the 
structure of the disc is rather similar to that obtained in the 
Roche lobe overflow simulation. The impact of the incoming 
stream seems therefore to have a rather small effect on the 
structure of the disc. One may note that in our simulations 
the width of the stream at the inner Lagrange point is always 
smaller than the height of the disc. This, combined with the 
fact that we use a mass ratio of unity, may explain this 
conclusion. 

(4) Concerning mass accretion rates, although we could 
not reach a definite conclusion due to the long time required 
to attain a steady state, we observe that the ratio between 
the mass accretion and mass injection rates is about 0.5 
to 0.6. This is however very dependent on the numerical 
viscosity and the inner boundary condition. 
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Figure 11. The particle location map on the orbital plane in the 
case of the Keplerian disc simulation with 7 = 1.2 at time T ~ 2.3 



4 CONCLUSIONS 

We have presented our results of 3D SPH calculations of ac- 
cretion discs in two different situations, Roche lobe overflow 
and the evolution of a Keplerian disc. The results of our 
simulations clearly show the following: 

(1) Even in the 7 = 1.2 and 7 = 1.1 cases, discs do exist. 
Their height, which depends only on the radial distance, 
is rather large but totally consistent with the hydrostatic 
balance equation. Compared with other authors, we use the 
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